#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use Getopt::Long;
use parse_bl2seq;
use db_parameters;
#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<'USAGE';

USAGE:

    get_aligned_sequences_hmr.pl
                            [--cdna_instead_of_protein]
                            [--help]
                            IDs.data


        This script takes a file of tab delimited gene_ids via
        STDIN or on the command line and prints the aligned clustal sequence
        pairs to STDOUT
USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $cdna_mode = undef;
GetOptions('help' => \$help,
		   'cdna_instead_of_protein' => \$cdna_mode,
		   );

die $usage if ($help);
die $usage unless @ARGV;





# connect to genome database
my $dbh = connect_to_panda();




my @species = ('human_mouse', 'human_rat', 'mouse_rat');
my @ids = ([0,1], [0,2], [1,2]);

my @gene_ids;
my $cnt_sequences = 0;
while (<>)
{
	chomp;
	my @genes = split /\t+/;
	next unless @genes == 3;
	$cnt_sequences++;
	push(@gene_ids, [@genes]);
}

print STDERR "$cnt_sequences\t three way alignments\n";
my $sql_cmd_orig = <<"PL/SQLCMD"; $sql_cmd_orig =~ s/\t\n/ /;
SELECT gene_id1, gene_id2, sequence1, sequence2, cdna1, cdna2, gaps1, gaps2, align_start, align_len
	FROM TABLE
	WHERE
		gene_id1 = 'GI1' AND
		gene_id2 = 'GI2';
PL/SQLCMD

use constant ENSG1 => 0;
use constant ENSG2 => 1;
use constant PROT1 => 2;
use constant PROT2 => 3;
use constant CDNA1 => 4;
use constant CDNA2 => 5;
use constant GAPS1 => 6;
use constant GAPS2 => 7;
use constant START => 8;
use constant LENTH => 9;


print "CLUSTAL W\n";
for my $gi(@gene_ids)
{

	for my $i(0..2)
	{
		my $table_name = 'orthologs.'.$species[$i]. '_sequences';
		my $sql_cmd = $sql_cmd_orig;
		$sql_cmd =~ s/TABLE/$table_name/;
#		print STDERR "$ids[$i][0]\t$ids[$i][1]\n";
		$sql_cmd =~ s/GI1/$gi->[$ids[$i][0]]/;
		$sql_cmd =~ s/GI2/$gi->[$ids[$i][1]]/;
#		print STDERR $sql_cmd, "\n";
		my $rows = $dbh->selectall_arrayref($sql_cmd)			or die DBI::errstr();
		#my $sth = $dbh->prepare( $sql_cmd) 					or die DBI::errstr();
		#$sth->execute()										or die DBI::errstr();
		#$sth->fetchall_arrayref([]);
		next unless $rows->[0];
		my @data = @{$rows->[0]};
		my ($seq1, $seq2);
		if ($cdna_mode)
		{
			$seq1 = substr(expand_cdna_with_alignment_str($data[CDNA1], "-", $data[GAPS1]),
							$data[START] * 3, $data[LENTH] * 3);
			$seq2 = substr(expand_cdna_with_alignment_str($data[CDNA2], "-", $data[GAPS2]),
							$data[START] * 3, $data[LENTH] * 3);
		}
		else
		{
			$seq1 = substr(expand_prot_with_alignment_str($data[PROT1], "-", $data[GAPS1]),
							$data[START], $data[LENTH]);
			$seq2 = substr(expand_prot_with_alignment_str($data[PROT2], "-", $data[GAPS2]),
							$data[START], $data[LENTH]);
		}
		print $data[ENSG1], "\t", $seq1, "\n",
			  $data[ENSG2], "\t", $seq2, "\n";
	}
}



$dbh->disconnect();

